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A treatment of direct simulation Monte Carlo method (DSMC) as a Markov pro- 
cess with a master equation is given and the corresponding master equation is de- 
rived. A hierarchy of equations for the reduced probability distributions is derived 
from the master equation. An equation similar to the Boltzmann equation for single 
particle probability distribution is derived using assumption of molecular chaos. It 
is shown that starting from an uncorrelated state, the system remains uncorrelated 
always in the limit N — > oo, where N is the number of particles. Simple applications 
of the formalism to direct simulation money games are given as examples to the for- 
malism. The formalism is applied to the direct simulation of homogenous gases. It is 
shown that appropriately normalized single particle probability distribution satisfies 
the Boltzmann equation for simple gases and Wang Chang-Uhlenbeck equation for 
a mixture of molecular gases. As a consequence of this development we derive Birds 
no time counter algorithm. We extend the analysis to the inhomogenous gases and 
define a new direct simulation algorithm for this case. We show that single parti- 
cle probability distribution satisfies the Boltzmann equation in our algorithm in the 
limit N — > oo, — ► 0, At — > where is the volume k th cell. We also show that 
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that our algorithm and Bird's algorithm approach each other in the limit Nk — ► oo 
where is the number of particles in the volume Vfc. 



I. INTRODUCTION 

Direct simulation Monte Carlo method (DSMC) 1] is a standard method to solve the 
Boltzmann equation numerically. In this method one divides space into cells of volume Vk 
(k = 1, 2, 3, ...) and takes a large number (JV) of simulated particles (10 3 — 10 6 ) to represent 
real gas molecules. The time evolution of the gas for a short time period At is calculated 
in two steps. In the first step some pairs of particles in the same cell are chosen randomly 
and are allowed to collide without changing their positions. A collision is allowed with 
a probability proportional to uT, where u is the relative velocity and E is the total cross 
section. In the second step all particles are propagated without collisions for a time At. 

The method is invented by Bird and Bird introduced the method based on physical argu- 
ments. A seminal paper of Bird [2] gives somewhat heuristic arguments to justify its use to 
solve the Boltzmann equation. One variant of the method was derived by Nanbu 3J] starting 
from the Boltzmann equation. Also it appears that essentially the same stochastic algo- 
rithms for a homogenous gas were invented independently by people interested in using them 



as a pedagogica 



distribution. 



4] 



tool to demonstrate evolution of a gas toward Maxwell-Boltzmann(MB) 



5] 



6]. In order to represent time evolution of the real gas such methods 



should converge to the true solution of the Boltzmann equation in the limit of N — > oo, 

q n 

Vk — > 0, At — > 0. Convergence proofs were given by Babovsky|7] and Babovsky and Illner[8| 
for Nanbu's method and by Wagner[9j for Bird's method. 

The cited convergence proofs are very formal and they appear to be written for math- 
ematicians. In this paper we give a simple derivation of Birds no time counter algorithm. 
We also show that, in DSMC, appropriately normalized single particle probability distri- 
bution satisfies Boltzmann equation for simple gases and Wang Chang-Uhlenbeck equation 
for molecular gases and their mixtures. The language of this development is familiar to the 
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physicist from the well known BBGKY hierarchy. 

In the next section we develop a general formalism for direct simulation. In order to 
demonstrate usefulness of the formalism we apply it to some simple money games. In the 
third section we apply the formalism to homogenous gases and show that, if appropriate 
collision kernels are chosen, the one particle probability distribution obeys the Boltzmann 
equation for simple gases and the Wang Chang- Uhlenbeck equation for molecular gases and 
their mixtures. In the fourth section we derive DSMC algorithm for inhomogeneous gases. 
Finally in the last section we give a summary and discussion. 

II. DIRECT SIMULATION AS A MARKOV PROCESS 
A. The Master Equation 

Assume that we have an assembly of things we call 'particles'. Particles can be real 
particles in a gas or humans or anything you can imagine. There are N particles in the 
assembly where N is a very large number. Each member of the assembly can be in any one 
of the 'states' where states are labeled by the parameter //. For a real gas /j, can be velocity 
vectors and for an assembly of people /i can be the money in their pocket on bank account. 
The \i can be discrete or continuous and it can stand for a collection of indices that can be 
both continuous and discrete. For the rest of this section we will treat /x as a continuous 
index. Integration over \i is actually integration over the continuous indices and summation 
over the discrete indices that /i stands for. 

We play a stochastic game with this assembly. We randomly pick pairs of particles 
and force them to 'collide'. A collision is an event that the particles change their states 
with a prescribed probability. Suppose we picked particles with states /ia and The 
probability that they will end up with state labels /ic and \in in the volume d^cd^D is 
T(/j,a, A*b; He, fJ>D)dfJ>cdfj,D where T(/ia, Hb', He, Hd) is the collision kernel. Collision kernel is 
assumed to be symmetric 

T(/i a ,Hb;Hc,Hd) = T(^ c ,^ d ;^ a ,Hb), (1) 
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T(fi A ,fi B ;fi C ,fiD) = T(fi B ,fi A ;fi D ,fi C ). (2) 

Also the probabilities are normalized 

J T(/i A , n B \ nc, Hd) d^c dfj, D = J T{ji A , fx B ; He, Hd) d^A dfi B = 1- (3) 

We define N-particle probability distribution /^(/^i, /jl 2 , Hn', n) such that 
/^(/-ti, fj,2, UN', n)dfj,idfj,2, dfJ,N is the probability of finding the particles 1,2, ...,7V 
in the d/iidfj,2, d/iN phase space volume after the n th collision. Since the particles are 
identical the f^ N \ni, fx 2 , HN',n) is assumed to be completely symmetric 

/ (JV) (^i,...,^-,...,^i, ...,n N ;n) = f^ifi!, fit, fij, fi N ;n). (4) 

We define reduced M-particle distribution as 

f i - M \n 1 ,...,H N ;n) = j / (iV) (A*i> ■■■■> Hn\ n) d^ M +i d{i M +2, dpL N - (5) 

We will denote f^ M \ni, Hm', n) (M = 1,2, ...,N) as f^ M \[i;n) shortly. As a convenient 
notation we also define f^\fi A , fi B ;n) as 

fij M) (VA,HB;n) = / (M) (/ii,...,/ii = fi A ,...,fi j= fi B ,...,fi M ;n), (6) 

where fa and fa; are replaced with fi A and fx B in f^ M \fa, fj, M ;n). Examples are 

fS\nA,HB;n) = f(fas,fa,faA,fa,-,fav;n) (7) 
fS\^A,HB;n) = f{fi l ,fi A ,fi 3 ,n B ,fi 5 ,...,fi N ;n) (8) 

We are ready to start now. The equation satisfied by the f( N \fa n) is given by 

1 N N 

f {N) (n;n+ 1) = — — — ttX!Z! / fip{HA,HB;n)T{fi A ,fi B ;fi i , f i j )dfi A dfi B . (9) 
N{N -1) l=lj ^J 

The meaning of this equation is clear. If the last pair we collided is i,j molecules, the 
probability of having fa, fa pairs at the end of collision is the probability of having initial 
states /jl a , fas (represented by fif\faA, fas', n)dfi A d^ B ) multiplied by the probability of ending 
with fa, fa (represented by T(/i a ,ii b ; fa, fa)). The sum over i,j and the factor 1/N(N — 1) 
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takes care of the fact that all pairs (respecting order of the molecules) are possible with the 
probability 1/N(N — 1). The state of the system after n+1 collisions depends on the state 
of system after n collisions and the direct simulation game is a Markov process actually. 
The eq.((9]) is the master equation for this stochastic process. 

In order to see clearly how this equation is derived let us multiply this with dfiid^-.-dfiN- 
The left hand side is 

p N \fi; n + l)dnidfi 2 ---d^N (10) 

and it is the probability of the system being in the phase space volume dfiid^.-.dfiN after 
the (n + l) th collision. On the right side we have 

1 N N 

N ( N J fij Va, HB;n)T(fi A , fJ-B] Hi, fij)dfi A dfiBd^idfi 2 -..dfi N . (11) 

(Here the integration is over \x A and lib only) In order to interpret this lets us look at % = 1 
and j = 2 term. It is the following term 

1 



f {N) (fiiA, Hb, fJ>3, A*4,--, fi N )dfiAdu JB dfi 3 du J4: ...dfi 



N 



N(N-l) 

x [T(fi A ,n B ;ni, yL 2 )diiidyL 2 ) (12) 

integrated over \x A , In this form the terms under the integration are product of three 
probabilities. 1/N(N — 1) is the probability of choosing i — l,j — 2 pair. The second 
parenthesis is the probability of finding the system in du J Adu JB du J3 du J4i ...du JN phase space 
volume before the collision. The last parenthesis is the probability of taking particles one 
and two from d^AdiiB to d\i\d\ii interval after the collision. When integrated over fj, A , 
Hb this term becomes the probability of arriving in dii\d\i 2 ---d\XN phase space volume after 
(n + l) th collision via a collision between particles one and two. If all such term are summed 
over i and j we find the probability of probability of arriving in d/j,id/j,2---d^N phase space 
volume after (n + l) th collision which is the same as eq. ljTO]) . 
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B. Asymptotic Behavior of the Master Equation 

Let us introduce a short notation for state variables: 

X = (xi,X2 ) ...,Xff) dX = dXidX2--.dX]y 

Y = (y 1 ,y 2 ,...,y N ) dY = dy 1 dy 2 ...dy N ■ (13) 
Z — (zi, Z2 t ---, Zn) dZ = dzidz 2 ...dzN 

Then the Master equation can be written in the form 

f(X; n + 1) = J P(X, Y)f(Y; n)dY. (14) 

The P(X, Y) has iV(iV — 1) terms and each one of the terms contains N — 2 delta functions. 
For example i — 1, j — 2 term reads as 

_ 1 , T(x 1 ,X2;yi,y 2 )8{x 3 - y 3 )...S(x N - y N ). (15) 

The general expression for P(X, Y) is 

-i N N I N \ 

p(x, y)= EE n^, *j\ v» n 8{*k-vk)\ (is) 

I ) i=i y fc^ij y 

The P(X, Y)dX is the probability that the system jumps from Y to dX phase space volume 
after a collision. As can be seen directly from eq. (|T6|) it is also symmetric: P(X,Y) = 
P(Y,X). As a probability density it satisfies the normalization condition 

J P(X,Y)dX = J P(X,Y)dY = 1. (17) 

We will need convolution of P(X, K) shortly. Let us define W(X,Y) as 

W(X,y) = J P(X,Z)P(Y,Z)dZ (18) 

It is easily seen that W(X,Y) is symmetric (W(X,Y) = W(Y,X)) and it also satisfies a 
normalization condition 

J W(X, Y)dX = J W(X, Y)dY = 1. (19) 
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Now we are ready to discuss asymptotic behavior or the master equation. Let us form 
Jf 2 (X;n + l)dX as 

J f 2 (X; n + l)dX = JdX^J P(X, Y)f(Y; n)dY^j (J P(X, Z)f(Z; n)dZ^j (20) 
= Jw(Y,Z)f(Y;n)f(Z;n)dYdZ (21) 

We can also write / f 2 (X; n)dX as 

J f 2 (X;n)dX = J W{Y,Z)f 2 {Y)dYdZ = J W{Y, Z)f 2 (Z)dYdZ (22) 

which follows from eq. (ll9p . Using these two relations we can write the following 

J f 2 (X; n + l)dX - J f 2 (X; n)dX = J W(Y, Z)f(Y; n)f(Z; n)dYdZ (23) 

W(Y,Z)f 2 (Y)dYdZ 



1 



2 
1 

~2 

The right side can be written as 



W(Y,Z)f 2 (Z)dYdZ 



J f(X- n + l)dX - J f{X- n)dX = ~J W{Y, Z) (f(Y; n) - f(Z; n)f dYdZ. (24) 

Since W(Y, Z) is always nonnegative the expression on the right is always negative or zero. 
This means / f 2 (X; n)dX decreases after each collision. The decrease stops when f(Y; n) — 
f(Z; n) — for all Y and Z and this means f(X; n) must be a constant. The equilibrium is 
reached when f(X;n) is microcanonical distribution. 

There is a final point to be discussed here. The above argument proves that the proba- 
bility density in the direct simulation always converges towards microcanonical distribution. 
If the phase space is divided in separate regions such that collisions cannot take the system 
from one region to another then the above argument must be modified. If Y and Z belong 
to different regions then W(Y, Z) = and f(Y; n) — f(Z; n) = is not required. But if Y 
and Z belong to the same region then W^Y, Z) ^ and f(Y; n) — f(Z\ n) = is required. 
This means that f(X;n) must be a constant in each region asymptotically but they can 
be different constants. For direct simulation of a gas total energy and total momentum are 
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conserved and the system stays on a constant total energy-total momentum shell. Asymp- 
totically the f(X; n) will be constant on each shell but they will be different constant for 
different shells. 



C. The hierarchy of Reduced probability distributions 



If we integrate the master equation over d/j,M+i, A*m+2, •••,/^jv we obtain the equation 
/^(^ + D= (iV -y_-f- 1) /'">(^) (25) 



2(N-M) JW 



ZTn ^ / (AU, Mb; T(/i A , /i B ; m, ^c) dfi A dfi B dp, G 



M(M-1 N M M 



7jEE/ f iff* (pa, Pb\ n) T(pa, Pb] Hi, A*i) <*AM <W 



The f( M \n;n + 1) depends on f( M+1 \/j,;n) and this represents a hierarchy of equations 



similar to the well-known BBGKY hierarchy 
The first equation in the hierarchy is 

/«(A*;n + l) = (l-2/iV)/ (1) (/i;n) (26) 
2 /■ 

+ iV 7 f (^a, HB;n)T(/j, A , /j, b ; /i c , h) d/J, A d/i B d/j, C - 

If we make the assumption of molecular chaos (AMC) 

f 2 \^ A) fx B ; n) = fM(ji A - n) f^(ji B - n), (27) 

we obtain a nonlinear equation for /^(//; n) similar to the Boltzmann equation. 

From now on we will suppress the superscript (1) in /^(/x;t) wherever it does not cause 
confusion. Using the relation 

f(/i, n) = J f(/i, n) f(nc, n) T(/i A , fi B ; He, H) d/i A d/j B djj c , (28) 

which follows from Eq.([n]) and the normalization of f(nc) an d imposing the assumption of 
molecular chaos we can write eq. f[2"6"|) as 

2 f 

f(p,n + 1) = f(jp\n) + — J [f, f] T([m a , Hb] He, h) dfi A dji B d^ c (29) 
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[f, f] = fil^A, n) f(/i B , n) - f(ii c , n) /(//, n) (30) 

A second simplification occurs for large N. The 2/N appearing in eq. (!29|) is a small 
number and we can take r = 2n/N as a continuous parameter which we call the collision 
time. Then At = 2/N and [f(^;n + 1) — f(fi;n)] /At can be taken as df(n,r)/dr. The 
eq. (!29l can be written in either of the following forms: 



8t 
df{fJt,T 



J [/, /] T(fj, A , fj, B ] fj, c , y) dfx A d^B dfx c - (31) 

-f(fi) + J f(fx A )f (^B)T(fXA,^B;i^c,^)dnAdfx B dnc- (32) 

We will call the first equation in the hierarchy 'the first equation' briefly for the rest of 
the paper. In latter parts of this paper we will call the integral on the right side of eq.( |3TT) 
'the collision integral'. From now on we will also suppress the collision time r in /(/x,r) 
wherever it is convenient. 



D. Justification of assumption of molecular chaos 

The only thing in this paper that is not fully rigorous is the assumption of molecular 
chaos. In order to have assumption of molecular chaos valid from the beginning we must 
start from an uncorrelated state 

/ {Ar) (/ii,yU 2 , ...,fi N ;n = 0) = hQix) h(p2)....h(/j, N ), (33) 

which is what is done in direct simulations mostly. The master equation eq.Q should be 
used to justify AMC. For finite N the AMC is not strictly valid and the AMC should get 
better and better as — ► oo. For M/N << 1 the eq. ( 1251) is written as 

f M \n-n + 1) = (1 - 2M/N) f M \^n) + 0(1/A 2 ) (34) 
2 M 



/ fiM+h^A^B\n) 
i=i J 

x T(ji A , hb] Hi, He) dfiA dfi B df^c 



where 0(1/N 2 ) are the terms of order 1/N 2 . If we invoke collision time r = 2n/N again 
and write [/ (M) (/i;n + 1) - / (M) (/!;n)j /At = df w {jiyr)/dr and we take the limit N -> oo 
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we obtain 




(35) 



dr 



c 



i=l J 

where M = 1,2,.. .,00. This is an infinite chain of coupled differential equations. If we 
invoke 



in the eq. fl35|) all the equations in the infinite chain are satisfies provided /W(/^;r) satisfies 
eq. (13T1) . This proves that in the limit N — > 00 the AMC remains valid for all r if we start 
from an uncorrelated initial state. 

What happens if we start from a correlated state that does not satisfy AMC? For finite iV 
there are always some correlations to any order. We know that the system evolves towards 
microcanonical distribution. In the limit iV — > 00 microcanonical distribution obeys AMC. 
This means even if we start from a correlated state the system will satisfy AMC better and 
better as the system evolves towards equilibrium for large N. Collisions destroys correlations 
and It should take only a few collisions per particle to destroy initial correlations. Moreover 
in the practical applications of DSMC in gas dynamics the N is almost always large and 
initial state is chosen as almost uncorrelated from the beginning. Therefore using the first 
equation to determine the single particle probability density is a justifiable process. 



We now show that expectation value (g(fi)) of a collision invariant g(/j.) is conserved. The 
g(fj) is a collision invariant if 



(36) 



E. Collision invariants and the H-theorem 



Ag = g(fj) + g{nc) - g{^A) - g{^B) = 0. 



(37) 



Multiplying eq. (j3Tj) and integrating over \i we obtain 



d 



/(//) g(fx) dfi= [f, f] g(n)T(n A , fi B ; fx c , n) dfi A dfj, B dfx c dji. 



(38) 
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Using symmetries of T(ha, Hb] He, h) an d relabeling integration variables among themselves 
we can write this as 

d If 

j- (9(h)) = 2 J [/> /] A # T (AU> Hb] He, h) d HA dfi B d Hc dfi. (39) 

The integral is zero because of eq. (157]) . 

We can derive an H-theorem for the first equation. Defining H(t) a 

H(r) = |/( M )ln(/(/x))d/x, (40) 

and using the eqs. f|T|2|) and eq. fl3Tj) we can express dff/ <ir as 
dH 1 /" 

~~d~T = ~\ \ T( ^> ^ ^) d ^ A d ^ B d ^ C d ^ ( 41 ) 

where 

m = If (Ha) I(hb) - f(H) f(Hc)} [Inf(HA) J(hb) - ln/(fi) /(/i C )] . (42) 
The $[/] can be shown to be always nonnegative as done in all kinetic theory books and 
T(ha, Hb', He, h) is intrinsically positive. Therefore dH/dr is nonpositive. There are two 
possibilities here. The H keeps decreasing toward negative infinity or it approaches an abso- 
lute minimum asymptotically and the system approaches toward an equilibrium distribution. 
Following the usual arguments of the H-theorem, the decrease of H stops only when 

ln/0u) + ln/( Ms ) = ln/( Mc ) + In /(//), (43) 

is satisfied which implies that In /(/i) is a collision invariant. If we choose the T{^a, Hb', He, H) 
such that there are collision invariants gi(fj) (i = 1, 2, L) then ln/(/z) must be expressible 
as a linear combinations of these collision invariants as 

ln f(H) = ci9i(h) + c 2 92(h) + ■■■ + c l9l(h), (44) 

where c%, cl are parameters describing the equilibrium. 

There is at least one trivial collision invariant. It is the number of particles entering and 
exiting the collision which corresponds to gi(n) = 1- When there are additional collision 
invariants the H has a lower bound usually. For the case of real gases momentum and energy 
are collision invariants and this makes H bounded from below. 
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F. Example: A game of discrete money gambling 

Here we give a simple example of a direct simulation money game with finite number of 
discrete states. Suppose everybody is given some random amount of money at the beginning. 
Everybody in the assembly has one, two or three dollars in their pocket. The random 
assignment of initial money ensures assumption of molecular chaos from the beginning. 
The collisions takes place as follows: Player 1 and player 2 share their total money such 
that nobody gets more than three dollars and both players get at least one dollar. All 
the possibilities satisfying these conditions have equal probabilities. If they have total two 
dollars (one dollar each) then the only possibility is that they will have one dollar each at 
the end with unity probability. If they have total three dollars then the possible outcomes 
are (1,2) and (2,1) with equal 1/2 probabilities. If they have total four dollars then possible 
outcomes are (1,3), (3,1), (2,2) with 1/3 probability each. If they have total five dollars then 
possible outcomes are (2,3) and (3,2) with 1/2 probability each. Finally if they have total 
six dollars (three dollars each) then the only possibility is (3,3) with unity probability. 

For this game the money is conserved in collisions and transitions between states with 
equal amount of total money is possible only. For N particles the total money can have 
values between iV to 3iV and there are a total of 2N + 1 separate regions in phase space. 
One cannot cross from one to another of these regions by making collisions. 

Now that we defined the game, how does single particle distribution evolves as we make 
collisions? The state variable /i is the amount of the money in the persons pocket and it 
takes the values 1,2,3. Let -P^(r) be the probability that a chosen person will have the money 
ji at the collision time r. From eq. (l32l) the P M (r) satisfies 

-P, + Pi T(l, 1, 1, 1) + P X P 2 T{1, 2; 2, 1) (45) 
+P 2 P 1 T(2, 1; 2, 1) + P 1 P 3 T(l, 3; 3, 1) 
+P3P1 T(3, 1; 3, 1) + P 2 P 2 T(2, 2; 2, 1), 



dPi 
~dr~ 



- = -P 2 + Pi T(2, 2, 2, 2) 



(46) 
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+P 1 P 2 T(l, 2; 1, 2) + P 2 P T(2, 1; 1, 2) 
+P 1 P 3 T(l, 3; 2, 2) + P 3 P T(3, 1; 2, 2) 
+P 2 P 3 T(2, 3; 3, 2) + P 3 P 2 T(3, 2; 3, 2), 



and 



-P 3 + P X P 3 T(l, 3; 1, 3) + P 3 Pi T(3, 1; 1, 3) (47) 



dPs 
dr 

+P3 2 T(3,3;3,3) + P 2 P 3 T(2,3;2,3) 
+P 3 P 2 T(3, 2; 2, 3) + P| T(2, 2, 1, 3). 

Inserting the T values this can be written as 

= "A + + ^ + f^iPa + \Pl (48) 

= -P 2 + -P| + PP 2 + -PxP 3 + P 2 P 3) (49) 

dP 2 1 

= -P 3 + -PP 3 + P 2 P 3 + P 3 2 + -Pi (50) 

This is a complicated set of nonlinear differential equations. But there are simplifying 
features because we know the collision invariants gi(fj) = 1 and g 2 (/i) = /x. Summing the 
eqs. (1481491150]) we obtain 

(Pi + P 2 + P 3 ) = (Pi + P 2 + P 3 - 1) (P + P 2 + P 3 ) , (51) 

dr 

and 

(Pi + 2P 2 + 3P 3 ) = (P + P 2 + P 3 - 1) (P + 2P 2 + 3P 3 ) . (52) 

The first equation tells us that since Pi + P 2 +P 3 = 1 at the beginning it always remains unity 
and probability is conserved. The second equation tells us that since Pi + P 2 + P 3 — 1 = 
always the expectation value (/jl) = Pi + 2P 2 + 3P 3 is conserved. 

We denote expected money in the pocket with m. We have two equations 

P1 + P2 + P3 = 1, (53) 
P + 2P 2 + 3P 3 = m, (54) 
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from which we solve P? and P^ as 



P 2 = _2P 1 + 3-m, (55) 
P 3 = Px + m-2. (56) 



Inserting P 2 and P3 in the eq. (1481) we obtain 



,/P| Pf + (m-^)P 1 + l(3-m) 2 . (57) 



rfr 1 v 3 ' 3 
Calculating roots of the quadratic term on the right we write this as 

dPi 



dr 



(Pi - n) (Pi - r 2 ) , (58) 



where r\ and r 2 are 



n = - ( 10 - 3/// V /M 3(/// - 11(3 - //;) ) . (59) 



r 2 = - \10 - 3m - yjl + 3(m - 1)(3 - m)j . (60) 

Notice that since 1 < m < 3 the term under the square root is always greater than or equal 
to unity. 

Solving eq.( l58i) is straightforward and we obtain 

D /x r 2(Po -r x ) -ri(p -r 2 )e- Ar 
(Po - ri) - (po - r 2 )e Ar 

where po = P\{t = 0) and A = r\ — r 2 . It is easy to verify that P\ (00) = r 2 and 
Pi (r) approaches this limit exponentially fast. One can check from eq. flBTJ]) that r 2 = 1 
at m = 1 and r 2 = at m = 3 and it behaves as it is expected. 

The conditions < P 2 < 1 and < P3 < 1 together with eqs. (l55|56l) gives conditions 
that .Pi(t) must satisfy. These conditions are expressed as 2 — m < P± < (3 — m)/2 when 
m < 2 and < Pi < (3 — m)/2 when m > 2. Therefore Pi(r = 0) initial value should obey 
these limitations. 

To find the equilibrium distribution directly without solving the differential equation we 
set dP^/dr = for fi — 1,2,3 in eqs. ()48|49|[50l and we obtain a set of algebraic nonlinear 
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equations. Setting P 1 = a, P 2 = ab, P 3 = ab 2 all three equations are satisfied provided the 
normalization condition 

a(l + b + b 2 ) = l, (62) 

holds. We were able to guess this solution from the H-theorem. There are two collision 
invariants g± — 1 and g 2 (fi) = A 4 - The second one is a result of conservation of money in the 
collisions. Therefore according to the H-theorem we must have lnP^ = C\ + C 2 fi and this 
gives the solution P^ = ab 11 ' 1 . We need one more relation to determine both a and b. This 
comes from expected money in the pocket: 

m = a (l + 2b + 3b 2 ) , (63) 

which is a conserved quantity during the 'time' evolution and it is set by the initial conditions. 
Solving these two equation we obtain 

a = - (lO - 3m - y/l + 3(m- l)(3-m)) , (64) 
b = (m-2 + v /l + 3(m-l)(3-m)^) /2(3 - m). 

Notice that a = r 2 and this agrees with solution of the differential equation. 
The H-function 

H = Pi In Pi + P 2 In P 2 + P 3 In P 3 , (65) 

is bounded from below for this problem since the function x In x is bounded from below and 
< P M < 1. We minimize H with the constraint that the expected money is fixed and 
probabilities are normalized. The constraints can be adopted with Lagrange multipliers. 
Taking the auxiliary function 

* = P 1 \nP 1 + P 2 \nP 2 + P 3 \nP 3 (66) 
-A 2 (P + P 2 + P 3 - 1) - A 2 (P + 2P 2 + 3P 3 - m), 

and setting d^/dPi = d^/dP 2 = <9\I>/<9P 3 = we obtain the same solution P M = ab^ 1 
where a and b satisfies the eqs.( !62ll63|) . The minimum value of H becomes 

H = a\na + ab\nab + ab 2 \nab 2 = ln^" 1 ' 1 ). (67) 
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G. Example2: A game of continuous money gambling 

Here we give another example of direct simulation money games with continuous states. 
In this case we were not even able to solve one particle probability distribution. We just find 
the equation for one particle distribution and guess the stationary one particle distribution 
from the H-theorem. We then show that it satisfies the equation for single particle probability 
equation. 

This time initially we give players a random amount of money between zero and, say, 
ten dollars. Suppose we pick a pair to collide, playerl has fii and player2 has \X2 amount 
of money. A computer produces a random number p between zero and one. Playerl takes 
p(/ii + H2) and player2 takes (1 — + /U 2 ) amounts of money and we pick another pair 
to collide. What is the final distribution when the system comes to equilibrium? 

The probability distribution that a person will have money /i satisfies the eq.( |32|) 

— - — = — /(/i) + / da db f(a) f(b) T(a, b, /i, v) da db dv, (68) 
or Jo Jo 

where the collision kernel is 

T(a, b, //, v) = —^—rS(a + b-/i- u)Q(a) Q(b) 9(/i) Q(u). (69) 
Here Q(x) is the standard step function 

x < 

0(x) = • (70) 

I 1 x > 

If we insert the T(a, b, /1, v) given in the eq. fl69l into the eq. fl68|) and perform the v integral 
we obtain 

«f = -/W + /\r*e ( » + t-^. (7i) 

ot Jo Jo a + 

This can be further simplified by changing variables x = a + b and y = a which yields 

a -M = . Htl)+ r dx r dy f(v)^-v) . (72) 

OT Jfj, Jo X 

The H-theorem insures that this equation will converge to an equilibrium distribution as 
t 00. Since we have money conservation in the collisions there are two collision invariants 



17 



gi(fj) = 1 and g^ii) = \i. Then the equilibrium distribution is 

fM = Ae- B ". (73) 

If the average money initially given to each person is m, the /(//) should satisfy two condi- 
tions 



o 



f(ji)dfji=l, (74) 

/ fj, f (/j) d/jL = m, (75) 
Jo 

and they fix the values of A and B in the eq. f)73p . The solution is 

fM = ^e~^ m . (76) 

ill 

If we insert this solution into eq. (|72|) we can easily check that right side of the equation 
becomes zero which confirms that f eq (n) is the equilibrium distribution. 

III. APPLICATION OF THE DIRECT SIMULATION FORMALISM TO 

HOMOGENOUS GASES 

A. Center of mass frame 

In the following sections we will need some results from studying the collision in the center 
of mass frame. Instead of deriving them for each case separately we derive the relevant results 
once for the most general case in this subsection and refer to formulae derived here as needed 
in the following subsections. In the rest of the paper bold letters denote vector quantities. 

Particles with states fi A = and \i B = v# and enter the collision and particles with 
states iiq = Vc and /jld — v exit the collision. We define the center of mass (CM) coordinates 
as 

H = (m A VA + rn B v B )/{m A + m B ) (77) 
H' = (m A w c + rn B v)/ (m A + m B ), (78) 
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and 

u = va — Vb, u = lul , n = u/w 

7 (79) 

u'= v c — v, u' = \u'\ , n' = u'/u' 
where is the mass of particles A and C and is the mass of particles B and -D. For 
one kind of gas all masses are equal and formulae for CM velocities H and H' reduce to 

H=(v4 + v B )/2, H' = (v c + v)/2. (80) 

Integrations over v^ and vg can be carried over in the variables H and u. The transformation 
between these two sets of variables are linear and the Jacobian is unity. Therefore 

/ /(vA,v B )rf 3 v^ 3 v B = J /(H,u)rf 3 Hrf 3 u. (81) 

In the following subsections we will deal with integrations over v A , v#, vp. Integrations 
over va, v b will be converted to integration over H and u in the CM frame. In each case 
there will be a Dirac delta function removing the integral over H. Integration over v c will 
be converted to integration over u' since v<^ = u'+v and there is no integration over v. 
Furthermore integrations over u' will be carried in spherical coordinates as 

J /(u'KV = J f(u')(u') 2 du'dn' (82) 

and in each case there will be a Dirac delta function removing the integration over u' . In 
the final expressions the integration over solid angle n' and u remain at the end. 

In order to evaluate the integrals we will encounter in the following subsections we must 
express v^, vb,vc in terms of the variables v, u, n'. This is a simple exercise in collision 
kinetics. We will do this for the inelastic collisions with unequal masses. This is the most 
general case we will deal in this paper. We will assume that molecules have internal energies 
e(A), e(B) and e(C), e(D). Let e = e(A) + e(B) and e' = e(C) + e(D). From energy 
conservation we have u'{u) = \Ju 2 + 2(e — e')/m r where m r = mArn B /{rriA + m B ) is the 
reduced mass and tua, are masses of the colliding particles. We can write u' = u'(u)n' 
and vc = v+u'(u)n'. From CM velocity conservation we have 

m A v A + m B w B = m A v c + m B w ={m A + m B )w + m A u'(u)ia! (83) 
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and we also have v A — v B = u. We solve v^, v B} v c from these as 



v A 



m A + m B 



-u'(u)n' + 



Vfi = v + 



-u (u)n 



m A + m B 
m A 



-u 



+ "IB 



-U 



u 



(84) 
(85) 
(86) 
(87) 



m A + m B 

vc = v+u'(u)ri 
l\u) = y/u 2 + 2(e- e')/m r 

For one kind of gas (m^ = m B = m ) without internal states (e(A) = e(B) = e(C) = 
e(-D) = 0) these equations reduce to 

v A = v + (un' + u)/2 (88) 

v B = v + (W-u)/2 (89) 

vc = v+wn' (90) 

Again for one kind of gas (m A = m B = m and m r = m/2) with internal states 
eqs. ()84ll85|l86|l87j) reduce to 

v A = v+ [u'(«)n' + u]/2 
v B = v+[u'(u)n'-u]/2 



vc = v+u'(u)n' 



u'(u) = \Ju 2 + 4(e - e')/m 



For a mixture of gases without internal states eqs. f)84ll85|86ll87l) reduce to 

m A 



v A = v + 



vb = v + 



/ . m B 
• un H u, 



"^A + m B m-A + TUB 

m A m A 
un 



Tn A + m B 



m A + m B 



u. 



vc = v+un . 

And for a mixture of gases with internal states eqs. f)84ll85ll86|87l) are the formule. 



(91) 
(92) 
(93) 
(94) 

(95) 
(96) 
(97) 



B. One kind of gas without internal degrees of freedom 

The state of particles are defined by three components of the velocity vector v. (We 
use bold letters for vectors throughout this paper) Bird's original algorithm to keep track 
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of time in the simulation was the 'time counter method'. Later Bird introduced 'No time 
counter method' (NTC) and declared time counter method 'obsolete' in his book. 1] Time 
counter method is more difficult (if not impossible) to formulate in the direct simulation 
formalism given in this paper and since NTC is the algorithm currently used we will derive 
NTC algorithms only in this paper. 

Here the state index /i refer the velocity vectors and the integration over /i stands for 
three integrations over components of velocities. The NTC kernel S(v A , v B ; v c , v) = Si + S 2 
is given by 

5 1 = ^5 (H - H') 5 (u 2 - (u') 2 ) a(n, n') (98) 

R ^ ' 

5 2 = (l - ^) S (v c ~ va) 5 (v - vfl) (99) 

Here cr(n, n') is the differential cross section and Sis the total cross section which is given 
by 

E = J a{n, n') dri, (100) 

where dn' is the solid angle in the direction of n'. The cr(n, n') depends on the angle 9 
between n and n' (n' ■ n = cos^). Hence <r(n, n') = a(n', n) and the kernel is obviously 
symmetric. The term 5(u 2 — (u') 2 ) = 5(u — u')/2u represents energy conservation and 
5 (H — H') represents conservation of center of mass (CM) velocity which is the same thing 
as the conservation of momentum. The kernel satisfies the normalization condition 

J S(v A , v B ; v c , v) d 3 v c d 3 v = J S(v A , v B ; v c , v) d 3 W d 3 W = 1. (101) 

Here the integral is taken in the CM coordinates. The Jacobian of the CM transformation 
is unity and d 3 u' = (u') 2 du'dn' . 

The S2 part of the kernel directly transfer initial velocities to the final velocities with a 
probability (1 — uL/ R) and hence causes a null collision. A null collision is a collision that 
particles do not change their states. The probability of making a real collision is 

J S 1 (v A ,v B ;v c ,v)d 3 v c d 3 v=^- (102) 
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where integral is calculated in the CM coordinates. 
Inserting 5(v^, v B ; v c , v) in eq. fl3TT) we obtain 

= /[/,/] SiKvBjvcvJdVAfltfvc (103) 

where 

[/,/] = /W/N-/(v c )/(v), (104) 
The S2 part of the kernel gives zero contribution in the collision integral 

/ [/, /] S (v c - v A ) S (v - v B ) rf 3 v A rf 3 v B d 3 v c = 0. (105) 

We evaluate the integral in eq. (11031) in the CM coordinates. We write d 3 v A d 3 v B = 
d 3 Hd 3 u and d 3 vc = d 3 u' = (u 1 ) 2 dv! dn'. When we do the integral we obtain 

df(y) 



^[/./lu^nVudn', (106) 



dr R 

where v^, vb, vp are expressed in terms of the variables v, u, n' in eqs. (l88|l89|90!) . 

The equationf ll06l) is essentially the Boltzmann equation with the difference that the 
Boltzmann equation is written for density in physical space. To obtain the Boltzmann 
equation we write this equation for F(y) = (N/V) /(v) where Vis the volume of the gas. 
Then we obtain 

dJ frT = R~G)J ^ F ^ F ^ ~ F (vc)F(v)} ua(n, n') d 3 u dn' (107) 

Now, if we change to the variable t = tV/RN = 2nV/RN 2 we obtain the Boltzmann 
equation for a homogenous gas 

= J [F(v A ) F(v B ) - F(v c ) F(v)) ua(n,n')d 3 udn' (108) 

Here t must be interpreted as the physical time and t = 2nV/RN 2 formula connects the 
physical time t and number of collision attempts n. 

Let us state the algorithm for a homogenous gas. We choose a number R big enough such 
that for only very few (say less than one in thousand) pairs uT*/R will exceed unity. We make 
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n = RNH/2V collision attempts to reach the desired time. For each pair we take a random 
number r and we allow the collision to happen if r < uTi/R. If the collision is allowed, we 
choose the direction of scattering n' according to the probability density <r(n, n')/E and a 
few more random numbers are used for that. Then we calculate and store final velocities for 
the colliding pairs and pick another pair. We keep taking and colliding pairs until we reach 
the desired time. 

Suppose the formula n = RN 2 t/2V yields 234.783 collisions. How do you make 0.783 
collisions? The way to do this in practise is to make 234 collisions first. Then throw a 
random number r and if r < 0.783 then go on to make a collision attempt. This can be 
justified from the formula 

2 f 

f{fi;n + 1) = f(n',n) + — J [f, /] T(fi A , fi B ; fi c , fi) dfi A dfi B dfi C - (109) 

After making n collision attempts with the NTC kernel S(v A , v B ; v c , v) we can change the 
kernel to 

P(va, Vb; Vc,v) = qS(v A ,v B ;v c ,v) + (I - q)5(v c -va) 5 (v - v B ) . (110) 

This kernel makes a NTC collision attempt with a probability q (which was 0.783 in the 
above example) and a null collision happens with the probability 1 — q. We use this kernel 
for the (n + \) th collision attempt (it is permissible to change the kernel) and this causes 
another At = 2q/N collision time and At = q(2V/RN 2 ) real time increase. 

C. Mixture of gases without internal degrees of freedom 

The state of particles are defined by three components of the velocity vector v and one 
kind index for which we will use p, q, r, s characters. We have M kind of gas without internal 
states in the mixture and there are N p number of molecules of the p th kind. The mass of p th 
kind molecule is m p . The probability density f(fi) = /(v,p) will be written as / p (v). 

Particles with states /j,a = (va,s), and /jl b = (v B ,r) enter the collision and particles 
with states fic = (vc,q) and ji B = (v,p) exits the collision. The integration over \i such as 
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/ f p { v )d^ stands for three integrations over v and summation over p. The center of mass 
(CM) coordinates are defined in eqs. (|77|78|79p . 

The NTC kernel G"(v^, vb; vc, v) = G\ + G2 is given by 

G l = |«J(H-H') 5(u 2 -{u') 2 ) a pq (n,n')5 pr 5 qS} (111) 
G 2 = (l-^) * (vo - v A ) * (v - v fl ) W (112) 

Here <7 p(? (n, n') is the differential cross section between gases of the p* h and q kind and 
E p(? is the total cross section which is given by 

S P9 = J a pq (n, n') tin', (113) 

where dn' is the solid angle in the direction of n'. The S pr 5 qs term in the kernel insures that 
particles do not loose their identities during the collisions. Again <r pq (n, n') = cr rs (n, n') due 
to the S pr 5 qs term and we also have the symmetry a pq (n,n') = a qp (n',n). The kernel is 
obviously symmetric. The term S(u 2 — {u') 2 ) and 5 (H — H') have the same meanings as 
before and the kernel satisfies the normalization condition 

M m . 

EE/ G r pq (v A ,v B ;v c ,v)d\ c d\= 1. (114) 
p =i q =i 

Again G2 part of the kernel directly transfer initial velocities to the final velocities with 
a probability 1 — (uYt rs )/R and hence causes a null collision. The probability of making a 
real collision is 

M M y 

EE (GiSKv B !V Cl v)<l 3 v c( l 3 v=f , (115) 

p=l q=l J K 

where integral is calculated in the CM coordinates. 

Inserting G r pq (yA, v B ; v c , v) in eq. (|3TI) and doing the summations over r, s and doing the 
integrals in the CM coordinates we obtain 
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where 

IF, F] = fK) F(y B ) - f%v c ) Fiy) (lis) 

Again we write this equation for F p {y) = (N/V) f p (v) and take t = 2nV/RN 2 to obtain 
Boltzmann equation for a mixture of homogenous gases without internal states 



dF p (v) 



M . 

E J l Fq M F p (v B ) - F«(v c ) F p (v)} ua pq (n, n') d 3 udn'. (119) 



Here v A , v B , vc are expressed in terms of the variables v, u, n' in eqs.( 

The algorithm is the same. We take n = RN 2 t/2V pairs and allow each collision with a 

probability (uS r8 )/R. If the collision is allowed we choose the scattering angle according to 

the cr rs (n, n')/£ rs probability distribution. 

Note that the normalization of / p (v) is given by 

M r 

E / / P (v)rf 3 v = l. (120) 
P =i J 

The integral / / p (v)<i 3 v is conserved during the simulation. From eq. (11161) its rate of change 
is 



J F(y) rf 3 v = / rf 3 v = E/ G p P l(v A , v B ; v c , v) (121) 

x [/ 9 (va) FM - f q (v c ) / p (v)] d 3 v A d 3 v B d 3 v c d 3 v. 



From normalization of probabilities in eqs. (|3llll4l) we have 



J G™(v A , v B ; v C) v) d\ c d\ = 1 (122) 
/ GH(v A ,y B ;v c ,v)d 3 v A d 3 v B = 1. (123) 



Using these relations the integral on the right side of eq. f)12ip can be written as 

M 

E 

9=1 



d 

dr 



J /P( v )d 3 v = ]T J p(v A ) F{v B )d 3 v A d 3 v B (124) 

9=1 

M r 

-Ey /Vc) / p (v)d 3 v c d 3 v. 
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These two terms are equal and they cancel each other yielding constancy of / / p (v)c/ 3 v. 
The number of molecules of the p th kind is 

N p = Nj f(v)d 3 v, (125) 
and it remains constant as it should. Hence the F p (v) is normalized as 

J F p (v)d 3 vd 3 x=N p} (126) 
where x is position of the molecule. 

D. One kind of gas with internal degrees of freedom 

For a homogeneous gas with internal states the \x stands for velocity v and a discrete 
index (for which we use a,/3,i,j) defining the internal quantum state of the molecule. The 
mass of the molecules is m. Particles with states ha = ( v a,/3) and hb = (ve,a) enter 
the collision and particles with states \ic = (yd) and fj,£> = (v,i) exits the collision. The 
integral over /i stands for integration over v and summation over the internal state index. 
The internal energy of molecule in the state 7 is E 1 and e = E a + Ep and e' = Ei + Ej. The 
center of mass (CM) coordinates are defined in eqs. (179f80l) . 

Let us define the no time counter (NTC) kernel K°f(\A, v b; v c ; v ) = Ki + K 2 where 

2u 



K 1 = h{n-H!)8 



and 



A e + M 2_ A e '_( u ')2 



K * = I 1 - 15 E E u ^tf ] 5 ( v c - va) <*(v - v fl ) S ia S jP . (128) 



u 



, <(n,nO, (127) 



Here m r = m/2 is the reduced mass where m is the mass of the molecules and R is a chosen 
parameter. The a°f(n, n') is differential and the E^- is the total cross section into the 
internal states i,j 

< = faff(n,n')da', (129) 
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where dn' is the solid angle in the direction of n'. This 
reciprocity relation of the inelastic scattering cross sections 



n 2 <(n,n') = (n') 2 ^(n',n), 



ternel is symmetric due to the 



(130) 



because {u/u') af? = (u'/u) a l h. 



The K 2 part of Kff{\A, Vs; Vc, v) directly transfers initial state to the final state and 



causes a null collision. The probability of making a real collision into the states is 

r u T, af3 

Pij = J Kt dv c dv =—£-. (131) 



R 



Therefore total probability of making a real collision is (J2iJ2j ^E^- )/R. 

Inserting the KfP (v^, v#; v<7, v) into the eq. (l3Tj) and doing the integrals in the CM 
coordinates we obtain 

9fliy) -^EEE / [//iN/ a N-/ j (v C )/ j ( V )]M?(n,n') ( i 3 udn'. (132) 



dr 



P 3 



"^ = EEE J [F p (y A ) F a {v B ) - Ffrc) F t (v)] ua°f(n, n') d 3 udn', 



(133) 



Here the part does not contribute to the collision integral as before. 

Again defining time as t = rV/RN = 2nV/RN 2 and defining the new functions i^(v) 
(N/V)fi(v) this is expressed as 
OF = 

a p j 

where Va, v b, vc are expressed in terms of the variables v, u, n' in eqs. fl91|92|93ll94p . These 
equations are the Wang Chang-Uhlenbeck equations for a gas with internal degrees of free- 
dom. Here the states are assumed nondegenerate for simplicity. Generalization to degenerate 
states is also very straightforward. 

Again we choose a number R big enough such that for only very few (say less than one 
in thousand) pairs u ^?j) / P will exceed unity. We chose n = RNH/2V random 

pairs. For each pair we take a random number r and we allow the collision to happen if r < 
J2j u ^tj K collision is allowed we choose the final state with the probability 



/{J2i J2j E£f ) and another random number is used to choose the final state. Finally we 



J 3 13 



choose the direction of scattering n' according to the probability density <r"- (n, n')/££- and 
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a few more random numbers are used for that. Then we calculate and store final velocities 
and state indices for the colliding pair and go on to choose next pair. 



E. Mixture of gases with internal degrees of freedom 

This case is a combination of previous two cases and it is very straightforward but unfor- 
tunately there are too many indices. The state of particles are defined by three components 
of the velocity vector v and one kind index for which we use p, q, r, s and one internal state 
index for which we use i,j, a, f3. We have M kind of gas with internal states in the mixture 
and there are N p number of molecules of the p th kind. The internal energy of i th internal 
state of p th kind molecule is Ef. The probability density /(//) = f(v,i,p) will be written as 
/f(v). 

Particles with states /i^ = (v^, f3, s), and /x# = ( v b, a, r ) enter the collision and particles 
with states nc = (y c ,j, q) and nr> = (v,i,p) exits the collision. We also define e = E% + E r a 
and e' = Ej + Ef. The integration over /x such as / ff (v)d/x stands for three integrations 
over v and summations over % and p. The center of mass (CM) coordinates are defined in 
eqs. (1771781751) . 

The NTC kernel is QfjffCvA, v b] v c ; v ) = Qi + Q2 where Qi and Q 2 are defined as 



Q, = 1<J(H - H') 6 
K 



2u 



jL e + u 2 - —e'-{u'f 



, oS5f(n,n')VV ( 134 ) 
u 



and 



Q2 = ( 1 - X) u ^iffgJ 5 ( v c _ v a) <^(v - v B ) 5 ia 8jp 5 pr 5 qs . (135) 

The delta functions 5 pr 5 qs insures that the molecules do no change identities during the 
collision. Here m r = mjjriB / {ttia + mg) is the reduced mass, R is a chosen parameter. The 
a°fpg q (n, n') is the differential cross section between species of the p th kind in the state a 
and q th kind in the state (3 and S^fpq 9 is the total cross section into the channel 

s Sf = / a Sr( n ' r/ ) dn ' ( i3e ) 



2s 



where dn' is the solid angle in the direction of n'. The Qtfp^i^A, v b; v cs v ) is also symmetric 



due to eq. (U30p . The Q2 directly transfers initial states to the final states and causes a null 
collision. The probability of making a real collision into the states is 

Pl3 = J Q ldwcdw= ^M^ (137) 

Therefore total probability of making a real collision is (J2%J2j u^lpq) I R- 

Inserting the Q^pqi^Ai v Bj v c> v ) into the eq. ([3"Tl) and doing the integrals in the CM 
coordinates we obtain 



dfliy) 



M 



EEEE / [A nf Q&tvA, v B ; v c , v) d^ A d\ B d 3 v c , (138) 

9=1 a & 3 

where 

[A P\t = flM - //(vc) ff (v). (139) 



After inserting Q^?~„ q we obtain 



A/ 



^§r^ = 4 E E E E / [A /1£f « ^(n, no ^ 3 u dn'. 



(140) 

The Q2 part does not contribute to the collision integral as before. Expressions of va, Vb, ve- 
in terms of v, u,n' are given in eqs. fl84ll85|86ll87p 

Again defining time as t = rV j RN = 2nV/ RN 2 and defining the new functions Ff (v) = 
(N/V) ff (v) this is expressed as 



dFF(v) 



M 

E E E E / F*{vb) - F«(v c ) i?(v)) (141) 



dt 

xu affpg Q (n, n') d 3 u dn'. 

These equations are the Wang Chang-Uhlenbeck equations for a mixture of gases with 
internal degrees of freedom. Here the states are assumed nondegenerate for simplicity again. 

Again we choose a number R big enough such that for only very few (say less than one 
in thousand) pairs (J2i J2j u ^Hj,pq)/R w iU exceed unity. We chose n = RNH/2V random 
pairs. For each pair we take a random number r and we allow the collision to happen 
if r < (J2iJ2juT l 'ff,p 1 g l )/R. If collision is allowed we choose the final state with the 
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probability S^^ 9 / (J2i J2j ^7j,pq 9 ) an d another random number is used to choose the final 
state. Finally we choose the direction of scattering n' according to the probability density 
a ij\pg 9 ( n ) n ')/^tf,pg 9 an d a few more random numbers are used for that. Then we calculate 
and store final velocities and state indices for the colliding pair and go on to choose next 
pair. 

Note that the normalization of ff (v) is given by 

EE//f(v)rf 3 v = l- (142) 

p i J 

The expression J2i I ff(y)d?v is conserved during the simulation. From eq. (jl38j) its rate of 
change is 

£e/ /fW 3 v = E/ ^d^ = j:Y,T,T,T, (143) 

Ul i J i J Ul q=l a p i j 

J [A f P }?f Q?f, P 7(v A , v B ; v c , v) d\ A d\ B d\ c d\ 

From symmetry and normalization of the kernel given in eqs. f|T|2|3|) we have 

EE/ QSVa, v b ; v c , v) d\ c d\ = 1 (144) 

» 3 

EE/ QtvV^A, v B ; v c , v) d\ A d\ B = 1 (145) 

Using this, we express eq.( ll43|) as 

d r m 

-r E / ft (v) rf 3 v = E E E / /^) /I(vb) d\ A d\ B (146) 

aT » J q=l a p J 

M . 

-EEE//;N/f(v)rf 3 v^ 3 v 

9=1 « 3 

These two terms are equal and they cancel each other yielding constancy of J2ifl{ v ) d 3 v. 
The number of molecules of the p th kind is 

N P = NJ2J /f(v)rf 3 v (147) 

i 

and as the above argument shows, it remains constant as it should. Hence the Ff(v) is 
normalized as 

E / F[(v) d\ d 3 x =N P , (148) 



30 



where x is position of the molecule. 



F. Relation to Kac's work 



Fifty years ago M. Kac 14| introduced a master equation similar to ours and derived the 
Boltzmann equation for a homogenous gas from it. Here we summarize his work and point 
out similarities. We will use a different notation than his. 

Suppose we have N particles in a gas contained in volume V. Collisions are as- 
sumed to take place randomly within the gas. Again we have a probability distribution 
/ (vi, v 2 , v^v; t) for their velocities. For brevity we will show this as /^(v; t) wherever 
convenient. Probability that the i th and j th particles having velocities and vg will collide 
and emerge with velocities vc-and vp in the phase space d 3 vc<i 3 V£> in a time interval dt is 
R(v A , v B ; vc, v D )d 3 vcd 3 v D dt. Here the R(y A , v b! v c> v d) is a function connected to differ- 
ential cross section but we will not need the precise relation until later. The total collision 
probability in dt time interval is S(y A , w B )dt where S^v^, v#) is given by 

S(v A ,v B ) = f R{v a ,v b ; v c ,v D )d 3 v c d 3 v D . (149) 

As usual we assume some symmetries for the R(y A , v B , v c , Yd) function: 

R(va, v b ; vc-, v d ) = R(-v c , v D ; v A , v B ), (150) 
R(v A , v B ; vc, v D ) = R{v B , v A ; v D ,v c ). (151) 

The /^(vi, v 2 , ...,Vjv;i) satisfies the master equation 



dfW(y) 



N N N N 



/ W (v)EE% v i) + EE / /g V) ( v A,v B )J2(v A ,v fl ;v i) v i )d 8 v >l d s v fl 

(152) 



^ *=i jVi »=i jV* 



In order to see where this comes from we write it for infinitesimal time interval dt: 

(N N \ 
i-*EE% v i) ( 153 ) 
<=i / 



/ N N \ 

+ dt \T,T, fiP(vA,VB)R(v A ,v B ;v i ,v j )d\ A d 3 v B 
\i=i j^i 7 / 
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Let us multiply both sides with d 3 v\...d 3 VN- Then f( N \v;t+dt)d 3 vi...d 3 V]y is the probability 
that the velocities are in the phase space volume d 3 v 1 ...d 3 VN at time t + dt. The first term 
on the right is 

(N N \ 
l-^EE^v,) . (154) 
i=l j# J 

The first parenthesis is the probability that the system was in <i 3 vi...<i 3 VAr phase space 
volume at time t and the second parenthesis is the probability that no collisions occurred in 
dt time interval. Their product is the probability of arriving d 3 vi...<i 3 V7v phase space volume 
at t + dt without making a collision. The second term in the right side are probabilities of 
arriving in d 3 vi...d 3 v^ by making collisions with different pairs. For example let us write 
i — 1, j — 2 term: 

/ (f {N) (v A ,v B ,v 3 ,...,v N )d 3 v A d 3 v B d 3 v 3 ...d 3 v N ) (R(v A , v b ; vi, v 2 )d\id 3 v 2 dt) . (155) 

The first parenthesis under the integral is the probability that the system was in the phase 
space volume d 3 v A d 3 VBd 3 V3-..d 3 VN at time t and the second parenthesis is the probability 
that the collision between particles one and two took them to (i 3 v 1 rf 3 v 2 phase space volume. 
If we integrate this product over v^, \ B we obtain probability of arriving in d 3 vi...d 3 \~N at 
time t + dt via a collision between particles one and two. To obtain total probability of 
arriving in <i 3 v 1 ...<i 3 vjv at time t + dt via a collision we sum such terms over all possible 
pairs. This argument clearly shows how the master equation is derived. 
Writing S^v^Vj) as 

S(vi,Vj) = J R(v A ,v B ;Vj,Vi)d 3 VAd 3 v B (156) 
the master equation can be written in a more symmetric form 



N N 



9t ^=l^ 



EE / (/fK,v B ) - /W(v)) A(va, v B ; v i5 ^)d 3 ^ A d\ B (157) 



All of the results we obtained from our master equation can be obtained for this master 

n 

equation too. Kac[14| showed that the distribution goes to microcanonical distribution as 
t — > oo. A hierarchy of reduced probability equations can be obtained for this master 
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equation too. Kac[l4( showed that in the limit ^ — oo if one starts from uncorrelated state 
at t — the system always remains uncorrelated. His arguments was different than ours. 
The first equation in the hierarchy (obtained by integrating over V2, V3, vjv ) is 

g/( ^ V) =2Nj (/ (2) (v a , v b ) - / (2) (v, v c )) R{v A , v B ; v c , w)d\ A d^ B d^ c (158) 

If we introduce AMC this equation becomes 

= 2Nj (f(y A )f(y B ) - /(v)/(v c )) R(v A , v B ; v c , v)d 3 v A d 3 v B d 3 v c . (159) 

Here the superscript (1) is dropped and time t is suppressed in /W(v;i). 

Now we go to center of mass frame (Equations 17911801) . In the CM coordinates the 
^( v A, v Bi v C) v ) is expressed as 

R{va, v B ; v c , v) = ^5 (H - H') 5 ( M 2 - (u') 2 ) a{n, n') (160) 

where V is the volume of the gas and cr(n, n') is the differential cross section. Inserting this 
into eq. (11591) and doing the integrals over the center of mass frame we obtain 
df(v) N 



dt V 



J [/(va)/(vs)-/(v c )/(v)] ua(n,n')d 3 udn', (161) 



where va,vb,vc are expressed in terms of the variables v, u, n' in eqs. ( |88ll89|90|) . If we 
write this equation for F(v) = (N/V)f(v) which is velocity distribution normalized to the 
number density per unit volume, we obtain the Boltzmann equation for a homogenous gas 

= J [F(y A ) F(vb) - F(y c ) F(v)] u a(n, n') d 3 u dn' . (162) 

Although both master equations have similar structures their philosophies are different. 
In Kac's work the collisions happens randomly and spontaneously in the gas whereas in direct 
simulation we take pairs and force them to collide. Direct simulation has applications to 
systems other than gases as we showed in the money games examples. In these systems there 
are not physical processes driving the collisions and instead we make the collisions. In Kac's 
work his motivation was to describe Boltzmann equation for g stochastic equation 

and the DSMC method had not been invented yet. Just as in our work, Kac's method can 
be generalized to molecular gases and gas mixtures and one can obtain Boltzmann equations 
for these cases by defining a suitable R(y A , v B ; vc, v) for each case. 
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IV. DIRECT SIMULATION FOR AN INHOMOGENEOUS GAS 

In this section we study NTC algorithm of DSMC method for inhomogeneous gas. We 
will not actually derive Bird's algorithm but we will define a similar algorithm to simulate 
inhomogeneous gas. We will show that single particle probability distribution of our algo- 
rithm satisfies the Boltzmann equation for an inhomogeneous gas. Then we will argue that 
both algorithms give the same results in the limit iV — > oo. 

We divide the physical space into cells as in the Bird's method. In our algorithm we 
take pairs not from the same cell but from all of the volume and we let each pair to make a 
collision attempt if both of them are in the same cell. 

We divide the physical space into cells and the k th cell has the volume V k . Now let us 
define the functions 

| 1 x G V k 

A fc (x) = { . (163) 

[ o x i v k 

We will also need the function 

r(x,x') = ]T kV tr ■ (164) 

k V k 

This function is zero when x and x' are not in the same cell and 1/V& when they are in the 
same cell. Its integral over x or x' is unity 

J T(x, x')d 3 x = J T(x, x')rf 3 x = 1. (165) 

At the end of this section we will take the limit 14 — ► 0. In this limit T(x, x') = for x^x' 
and T(x, x') = oo for x = x' and eq. (11651) is still satisfied. These are properties of the Dirac 
delta function and we have the limit 

hinr(x, x') = <J(x - x) (166) 

Now we can start the discussion. We will treat the simplest case for clarity. We develop 
our arguments for one kind of gas without internal degrees of freedom. The generalization 
to the other cases is very straightforward and will be summarized at the end of the section. 
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The state index /x represent position of the particle x and the velocity v. The collision 
kernel is Z — Z\ + Z 2 where Zi and Z 2 are 

Zi(x a v a ,x b Vb;x c Vc,x d v ) = S(v A , v B ; v c , v D ) (167) 

xT(x A , x B ) fi5(x c - x A ) 5(x D - x B ) , 

and 

^(x^jXbVbJXcVcXdVb) = (1 - fir(xA,x B )) 5(x c - xa) (168) 

x5(x £) - x B ) 5(v c - va) 8(v d - V B ). 

Here 

-1 

(169) 

is a constant chosen to insure that probability of making a collision in any cell is less than 
unity. The S(va, v b ; vc, v d ) is given in eqs. (I98|99j) . The Z 2 does not change states of the of 
the particles and the pair will not be allowed to make a collision attempt with a probability 
(1 — Q r(xA, x B )) . The probability of a collision attempt is fir(x A , x B ) and in a real collision 
positions of particles do not change because of the 5(xc — X4) — x B ) term in the Z. 
The Z(x j 4V j 4, x^v B ; xc-vc, X£)Vb) is symmetric and satisfies the normalization condition 

J Z(x j 4V J 4,x jB v B ;xcVc,x D v D )d 3 v j 4(i 3 v j B(i 3 x j 4(i 3 x j B = 1, (170) 
J Z(x A v A ,x B v B ;x c Vc,x D v D )d 3 Vcd 3 v D d 3 x c d 3 x D = 1. (171) 

The f2r(x j 4,x B ) vanishes unless x A and x B are in the same cell and f2r(xA,x B ) = fl/Vk 
when x A and x B are in the cell The probability of having both particles in the cell Vk is 
(Nh/N) 2 where Nk is the number of particles in the cell Vk during the collisions part of the 
simulation. Therefore the probability of a pair making a collision attempt in the k th cell is 
Pk = (fl/Vk)(Nk/N) 2 . The l/V^term looks awkward in this probability but it is absolutely 
necessary as the following argument shows. Suppose the physical density is uniform and 
therefore Nk/N = Vk/V where V is the total volume. When density is uniform we expect 
that the probability of having a collision in Vk is proportional to When Nk/N = Vk/V 
in inserted in Pfcwe find pk = QVk/V 2 which is proportional to Vk as expected. 
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Now we insert the kernel Z in the eq. (l3"Tj) to obtain 
df{xv, t) 



dr 



J [f, f] Z(x A v A , x B v B ; x c v c , xv) (172) 
xd 3 w A d 3 v B d 3 v c (fx a d 3 *B d 3 *c 



where [/, /] is 

[/,/] = /(xav a , r)/(x B v B , r) - /(x c v c , r)/(xv, r). (173) 

The Z 2 part of the collision kernel does not contribute to the collision integral. After doing 
the delta function integrals over positions x^ , x# we obtain 

df{ J r ' T) = O / [/(x'v,, r)/(xv s , r) - /(x'v c , r)/(xv, r)] (174) 
xr(x, x) S(v A , v B ; v c , v)d 3 v A d 3 v B d 3 v c d 3 x'. 

Now we insert 5 = S± + S% from eqs. (I98|I99I) in this equation. The S 2 part gives no contri- 
bution to the integral as before. Doing the integrals over v A , vb , in the center of mass 
coordinates we obtain 

dfi ™ ,T) = § / If&VA, r)/(xv fl , r) - /(xV C) r)/(xv, r)] (175) 
xT(x, x) o-(n, n) d 3 u dri d 3 x'. 

where v A , vg, vp are given in eqs. (l88f89ll9~Uj) . In order to have complete correspondence with 
the Boltzmann equation we define the new function F(xv, r) = iV/(xv, r) and we also 
define the new variable t = Qt/RN = 2Qn/RN 2 to obtain 

= LcF(xv, t ) (176) 

where the operator L c is defined as 

L c F{xv,t) = J [F{x!v A ,t)F{xv B ,t) - F{x?v c ,t)F{xv,t)] (177) 
x T(x, x') cr(n, n') d 3 u da! d 3 x'. 

Here t is interpreted as the physical time. 

In the collisions part of the DSMC method we make collision attempts for a time At where 
At is a small time interval. This corresponds to At = RNAt/Q collision time passage or 
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An = RN 2 At/2Q pairs chosen. From eg. (11761) . after making An collisions attempt F(xv, t) 
becomes F*(xv,t) 

F*(xv, t) =(1 + AtL c )F(xv, t)+0((At) 2 ) (178) 

where 0((At) 2 ) is an error term of order (At) 2 . 

Next we perform free propagation step where x — > x+Atv and v — > v+Ata transforma- 
tion is made for each particle. Here a = F/m is the acceleration of the particle due to the 
force F and it can depend on both position and velocity of the particle. This changes the 
N particle distribution function f( N \x.i, Vi; x 2 , v 2 ; ...,xjv, Vjv) to 

/ (JV) (xi - Atvi, Vl - Atai; x N - Atv N , v N - Ata N ). (179) 

The jacobian of the transformation is unity with a correction of order (At) 2 and there- 
fore this expression is correct with an error of the same order. Integrating this over 
x 2 , v 2 ; ...,Xjy, Vjv we find that the single particle probability distribution /^(x, v) changes 
to /W (x-Atv, v-Ata) with an error term of order (At) 2 . Therefore F*(x, v,t) becomes 
F*(x— Atv, v— Ata,t) which is taken as F(x, v, t + At). Hence 

F(x, v, t + At) = F*(x-Atv, v-Ata,t). (180) 

Using eq. (II 781) and expanding F(x— Atv, v— Ata,t) up to first order terms in At we obtain 

F(x, v, t + At) = ( 1 - Atv^- - Ata^- + AtL c ) F(x, v, t) + 0((At) 2 ) (181) 
\ dx dv J 

where 0((At) 2 ) is the error terms of order (At) 2 . Taking the limit At — > we obtain 

OFfr v,«) +v , aF(x,v, i ) + = ~ 

at ox m ov 

This equation is similar to the Boltzmann equation but it is not the same. Already 
when treating r = 2n/N as a continuous parameter we took iV — > oo limit implicitly. The 
remaining limit is V& — ► and we know that T(x, x') — > 5(x — x') in this limit. After setting 
T(x, x') = <5(x — x') performing the x' integral the operator Lc reduces to 

L c F(xv) = J [F(xv A , t)F(xv B , t) - F(xv c , t)F(xv, t)] a(n, n') d 3 u dn' . (183) 
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With this form of the Lc the eq. (1182j) is the Boltzmann equation. 

Hence we have shown that in direct simulation algorithm for inhomogeneous gas the one 
particle probability distribution satisfies the Boltzmann equation. Now, how do we connect 
this to the Bird's NTC algorithm? Clearly they are not the same. In fact our algorithm is 
not practical since great majority of chosen pairs will not be in the same cell and therefore 
will not make collisions. 

In the time interval At we choose An = RN 2 At/2Q pairs. The probability that each 
pair will make a collision attempt in the k th cell is p k = (Q/V k )(N k /N) 2 . Let n k be the 
number of collision attempts that take place in V k . The expected value of n k is 

RN 2 

n k = An-p k = -r^At. (184) 

This is the same as number of collision attempts in V k in Birds algorithm. The difference 
is that in Birds algorithm the number of collision attempts in each cell is fixed as n k = 
RN k 2 At/2V k whereas in our algorithm the n k has a probability distribution with a mean 
value RNlAt/2V k . The probability distribution for n k is given as 

(AnV 

P{n k ) = { ; v ( Pk T k (l-p k r n . (185) 
(An - n k )\ [n k )\ 

In the limit of V k — > we have p k — > and the P(n k ) becomes the Poisson probability 
distribution 

P(n k ) = ^^-exp(-n k ). (186) 

The width of distributions in eqs. (H85|ll86l) is of order ^JW^. For large values of n k we have 
n k /n k = 1 + 0(l/y/Wk) where 0{1/ yJWk) is a term of order l/y^. 

Now we take the limit N k — > 00 and 0{l/\/Wk) error term vanishes. In a more mathe- 
matical language, probability that n k /n k = 1 is unity. Hence both methods approach each 
other in the limit Nk — * 00 and single particle probability distribution in Bird's method too 
should satisfy the Boltzmann equation (eq. (11821) ) in this limit. 

There is an important distinction in the limits taken for both method to satisfy the 
Boltzmann equation. In our algorithm we take iV — > 00, At — > and V k — *■ limits. This 
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does not mean that number of particles in each cell (N^) will go to infinity. For example for 
a uniform density we have _/V& = (N/V)Vk- Here N oo and — > limits does not imply 
anything about Nk- NVj~ can remain finite and even can go to zero and still our algorithm 
satisfies the Boltzmann equation. The Bird's algorithm requires Nk — > oo to satisfy the 
Boltzmann equation however and this is a more stringent requirement. 

We did this analysis for the simplest case of one kind of gas without internal degrees of 
freedom for clarity. It is very simple to generalize this to the other cases by replacing the 
kernel S in eg. (1X671) with G™ in eg.(TTTTl) or with K°f in eg. (11271) or with Qff^ in eg.(TT34j). 
Then the Boltzmann eguation will be replaced by the Wang Chang-Uhlenbeck eguation but 
all of the arguments will remain the same. 

V. CONCLUSIONS 

Let us list our contributions in this paper. 

• In this paper we introduced a general formalism for direct simulation processes. We 
defined the direct simulation as a markov process with a master eguation and we found 
the master eguation given in eg.(??). Definition the DSMC algorithm as a stochastic 
process governed by a master eguation does not exist in the literature of the DSMC 
method to our knowledge. 

• Starting from the master eguation we showed that the N-particle probability density 
evolves towards microcanonical distribution as the number of collisions go to infinity. 

• We derived a hierarchy of eguations similar to the BBGKY hierarchy for the reduced 
probability densities given in eg. (l25l) 

• We showed that if AMC approximation is employed the single particle probability 
distribution satisfies an eguation given in eg. (l29|) . In the limit iV — ► oo this reduces 
to eg. (13 II) which is an eguation similar to the Boltzmann eguation. 
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We found the equations of the hierarchy in the limit N — > oo (the eg. (135]) )and showed 
that the ansatz /^(/xi, 1^2, •••> A*m! t ) — /^(/^li 7 ") f^Kl- l 2i T )----f (a*m;t) satisfies all 
the equations in the hierarchy provided the /^(/j;t) satisfies the eq.( l3TT) . This ensures 
that in the limit N —>■ 00 the AMC is satisfied for all times if one starts from an 
uncorrelated initial state. 

We gave two simple examples from direct simulation money games. The discrete 
money game example has the nice feature that it is exactly solvable and we observe 
from the solution that the approach to the equilibrium is exponentially fast. 

We obtained the H-theorem and conservation of expectation values of collision invari- 
ants. These results are familiar to most readers from the standard treatments of the 
Boltzmann equation. But it is worth repeating them here because although the equa- 
tions are similar they are applied to wide variety of different problems in the direct 
simulation setting, not just to gases. 

We applied the formalism to the direct simulation Monte Carlo method for real ho- 
mogenous gases which is a standard method to solve the Boltzmann equation. Intro- 
ducing appropriate kernels we obtained NTC algorithm for a homogenous gas and we 
showed that the appropriately normalized single particle probability distribution sat- 
isfies Boltzmann equation for simple homogenous gases and Wang Chang-Uhlenbeck 
equations for homogenous molecular gases and their mixtures. The derivation of con- 
servation of J / p (v) d 3 v for mixture of gases without internal degrees of freedom and 
J2i I ff ( v ) for mixture of gases with internal degrees of freedom should be also 
familiar to the reader from the standard treatments of the Boltzmann equation. The 
novel feature of our derivation is the significant simplification that the normalization 
of T(ha, Hb, He, A*d) given in the equations fl3"1ll22lll23|ll44lll45l) provide to obtain the 
result. If we try to obtain the same result from the Boltzmann equation we would 
have to use the argument that the integrals in (I122lll23lll44|ll45p are functions of the 
collision invariants. 
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• We introduced a new algorithm to do the DSMC calculations for an inhomogeneous 
gas. Our algorithm is not practical for the actual practice of the art because of wasting 
the great majority of the chosen pairs. We showed that the single particle probability 
distribution satisfies the Boltzmann equation in our algorithm in the limits N — ► oo, 
At — > and V& — > 0. We also showed that Bird's algorithm for DSMC converges to 
our algorithm if — » oo is taken in addition to the limits At — > and V% — > 0. Birds 
algorithm requires more stringent requirements to satisfy the Boltzmann equation. 
To prevent any misunderstanding we stress here that our algorithm is not intended 
as a practical scheme to implement DSMC calculations. The Bird's algorithm does 
not easily fit in the direct simulation formalism presented in this paper whereas the 
algorithm we presented does. We showed that our algorithm gives the Boltzmann 
equation in the limits iV — > oo, At —>■ and — ► and we also showed that our 
algorithm and Bird's algorithm converges to each other if we go to more stringent limit 
of Nk — > oo. Therefore we proved indirectly that Birds algorithm satisfies Boltzmann 
equation in the limit — > oo, At — > and — > 0. Therefore we introduced our 
algorithm as a tool to study convergence of Bird's method and not as a practical way 
of doing DSMC calculations. 

Meaning of the convergence here should be interpreted according to the ensemble theory 
of statistical mechanics. We imagine practically infinite number of identical systems (com- 
puters with human operators) doing the same direct simulation and call this the ensemble. 
The r)d/i represents ratio of number of particles in d\i to the total number of particles 

averaged over all the ensemble. When you perform a direct simulation on a computer you 
are just one member of the ensemble. Your results will show statistical fluctuations. But 
when you do the same simulation many times with different initial states chosen according 
to a uncorrelated probability distribution f^ N \fJ,±, jU 2 , I^n'i n = 0) = h(fi 1 )h(fj, 2 )----h(fiN) 
you form your own ensemble and averages over them will nicely follow f^(fi;r) obtained 
by solving eq. (l3il with the initial value /^(//; r = 0) = h(jj). 
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This work can generalize to chemical reactions and radiative processes in a more or less 
straightforward fashion. But there are enough number of subtleties such that we leave them 
to future publications. 

A simplified version of this paper [2] containing only one kind of homogenous gas without 
internal degrees of freedom is published in American Journal of Physics. The material in 
that paper makes a small fraction of the material in this paper. The present paper contains 
much new material and overlap between the two papers is small. 
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